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I. INTRODUCTION 

The Bethe-Peieiis or cavity method has a long his- 
tory in statistical mechanics pQ. The application of this 
method to disordered systems has recently undergone a 
considerable revival, mainly in connection with the anal- 
ysis of typical case complexity of random NP-complete 
(i.e. difficult) optimization problems. This recent work 
has led to an improved understanding of the statisti- 
cal mechanics of disordered systems in particular to a 
formulation of the physics of replica symmetry breaking 
without resorting to replicas. Importantly, it has also led 
to a new class of algorithms, now known as survey prop- 
agation, for the optimization problems In addition, 
very recent work [3j [4j (5] has attempted to generalize the 
cavity method/belief propagation to disordered quantum 
systems, obtaining encouraging results. 

The work alluded to above does not address one 
striking feature of the physics of disordered systems, 
namely the presence of Griffiths-McCoy (GM) singulari- 
ties O [TJ [8] in their thermodynamics in an applied field. 
While the form of these singularities can be readily deter- 
mined from rough estimates of the statistics of the rare 
regions from which they emanate, their detailed extrac- 
tion can be a tricky task due to their delicate nature, 
especially in classical systems. Indeed, in field theoretic 
formulations they appear as non-perturbative (instanton) 
effects 0. 

In this paper we consider the task of extracting these 
GM singularities from the cavity method. The method 
is exact on Bethe lattices and hence the functional re- 
cursion relation to which it gives rise must contain GM 
physics which exists already on these lattices. The chal- 
lenge then, is to extract it by constructing the appropri- 
ate fixed point solution. To this end we study the par- 
ticular case of a diluted ferromagnet on the Bethe lattice 
which has an extended GM region at low temperatures 
and large dilution. While the general, exact, determina- 
tion of GM singularities everywhere in the phase diagram 
is a hard problem, we are able to solve the problem in the 
infinite coupling limit made precise below. Here we can 



directly solve the cavity equations in a field and relate the 
solution to the statistics of clusters and to the density of 
Lee- Yang zeroes commonly used to characterize GM ef- 
fects. Further, in this limit the phase transition between 
the paramagnet and the ferromagnet is itself essentially 
of a GM character and its critical behavior, which we 
extract, can be viewed as an enhanced GM phenomenon. 

The problem of the dilute Bethe lattice ferromagnet 
and of GM singularities has been considered before us 
[TUl ITT] by different methods. As our interest is primar- 
ily in the development of the cavity method, we give a 
self-contained presentation in this paper from that view- 
point. We turn now to a more detailed enumeration of 
the contents of this paper. 

II. MODEL AND ORGANIZATION 

We consider the following disordered Ising Hamiltonian 
on the Bethe lattice with connectivity q: 

/3H e = -J^eijCTiO-j - H^cri, (1) 

(ij) » 

where 

_ J 1 with probability p; 
iJ 1 with probability 1 — p. 

The random couplings indicate the presence or ab- 
sence of a bond in the diluted Bethe lattice. For probabil- 
ity p < p c = l/(q—l) the lattice has no giant clusters and 
the density of large finite clusters decays exponentially. 
For p > p c , giant clusters exist with finite density and 
at the percolation transition, p = p c , the density of clus- 
ters of size n develops a long algebraic tail W n ~ n -5 / 2 
(independent of q) [T^]. We note that the dimensionless 
coupling constants J and H differ from the conventional 
magnetic exchange and field by factors of inverse tem- 
perature 0=1/T. The limit T -> with J > H will be 
denoted as the J = oo limit. 

In Section |III| we provide a guided tour of the well- 
known phase diagram of this model from the point of 



2 



view of the cavity method. We then establish the critical 
behavior at the phase transitions using a set of recursion 
relations for the moments of the cavity field distribution. 
We also show that the critical behavior can be extracted 
via a simple numerical algorithm, which we discuss in 
some detail in the Appendix. 

The following sections are devoted to investigating ex- 
act analytic results in the infinite (dimensionless) spin- 
spin coupling limit, J = oo. This corresponds to the hor- 
izontal axis of the phase diagram in Figure [l] In Section 
IV[ we find an explicit expression for the magnetization 
M(H) by means of a sum over connected clusters, which 
follows the standard GM treatment due to [T3]. In Sec- 
tion [V] we show that the same expression can, in fact, 
be extracted from the cavity method. In this limit, the 
magnetization goes to zero with the field for p < p c — 
l/(<7 — 1) while for p > p c a spontaneous magnetization 
develops. We first show that: 1) for p < p c the asymp- 
totic series expansion for the magnetization contains only 
integer powers M (H) = xH+C3H 3 +... and 2) on the con- 
trary at p = p c the series expansion contains semi-integer 
powers as well M(H) = c 1/2 y/H + ciH + c 3/2 H 3 / 2 + .... 
That is, the critical exponent 6 = 2 at p = p c , J — oo. 

In Section |VI| we develop an alternative integral rep- 
resentation for M(H) that corresponds to a harmonic 
expansion. This representation will allow us to calculate 
the (smoothed) density of Lee- Yang (LY) zeros p sm at 



J = oo on the imaginary H axis (8 — ImH) in Section 
VII and to show that for p < p c a GM phenomenon in- 
deed occurs, i.e. the density of zeros is non-zero and van- 
ishes as e~ a l® when approaching the origin. For p = p c 
we find a = and the density vanishes as the power law 
p cx \/6. 

Finally, the promised Appendix briefly describes the 
"population dynamics" algorithm used in the numerical 
work. 



III. PHASE DIAGRAM AND CAVITY 
EQUATIONS 

The p — J phase diagram (Figure [TJ of the diluted fer- 
romagnet Tl e is physically well understood and can be 
derived naturally in a cavity method formalism |14| 1X5] . 
In this approach, one considers the flow of cavity fields 
from the boundaries of the tree inward toward the cen- 
ter. A cavity field hi on a spin at a distance d from 
the boundary describes the spin's magnetization in the 
absence of the link connecting it to the next spin inward. 
The cavity field hi only depends on the cavity fields on 
CTi's neighbors at distance d—l and therefore one can de- 
fine a natural flow for the depth dependent distribution 
of fields P^(h): 



(h) = E e J f J] dh^-V (hi^j s(h - <K + H, Je t ) 



where 



(2) 



u(hi + H, Jti) — tanh 1 (tanh(Jej) tanh(fti + H)) 



(3) 



gives the bias on the field h due to a spin Ui connected 
through a link Jti. E e is the expectation with respect to 
the 6ij distribution. Fixed point distributions P^°°^(/i) 
describe the statistical features of the bulk (central re- 
gion) of the Bethe lattice. In order to break the Ising 
symmetry, we will always assume an infinitesimal uni- 
form positive boundary field P(°'(ti) = 5(h — + ) as the 
starting point for the flow. 



In the undiluted model, p — 1, all of this discussion 
reduces to the simple Bcthe-Peierls mean field theory for 
a connectivity q lattice. Since there is no randomness, 
the cavity field distributions P^ are simply delta func- 
tions located at, possibly depth dependent, fields h^ d \ 



Equation (|2j) reduces to a flow equation for h^: 

9-1 

h(d) = XX ft(d_1) + H ' J ) ( 4 ) 

i=l 

= (q-1) tanh -1 (tanh(J) tanh(/i (d-1) + H) ' 



For J < J G = tanh 1 f^zfji th° now at if = has 

only one fixed point = corresponding to the 

paramagnetic phase. For J > Jq, the ft/ 00 ) = fixed 
point becomes unstable to a spontaneously magnetized 
ferromagnetic fixed point with ft/ 00 - 1 > 0. Expansion 
of the fixed point equation to leading order in H and 
e = ( J — Jq) gives the well-known mean- field critical 
exponents at J = Jq: 



M(H,J=J G ) ~ JjV*; 



6=3 



M{H = 0,J>J G ) ~ {J-JgY; 13 = 1/2 



(5) 
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FIG. 1: Phase diagram for the diluted ferromagnet on the 
connectivity q — 3 Bethe lattice. 



Under dilution, we must return to the more general 
cavity distribution flow defined by equation ([2| to extract 
the phase behavior. Notice that the paramagnetic cavity 
distribution P PM (h) = 6(h) is always a fixed point of 
the flow at H = 0, just like = is always a solution 
for the undiluted model. As in undiluted case, this fixed 
point will become unstable above some critical coupling 
J c (p)- Near P PM (h) (i.e. for small h), we consider the 
linear stability of the first moment of P^ d \h): 



(h) 



E, 



E, 



J dh (jj dhiP^ihi^J 5^h-J2 u ^j h 

^ \i=i / i=i 

« ( g -l)(E e tanh(Je))(/ l ) ( ^ 1) 

to leading order. Thus, 1 = (q — l)ptanh(J c (p)) gives the 
critical boundary separating a stable paramagnetic phase 
from the ferromagnetic phase. A small rearrangement 
gives: 



J c (p) = tanh 1 ( — 



(6) 



This agrees precisely with the undiluted critical point 
J c (p — 1) = J a found above and also predicts that for 
p < p c the paramagnetic phase persists for all finite J. 
There is no ferromagnetic phase transition for a model 
with only finite clusters, as one expects. 

In order to extract the critical behavior along the di- 
luted phase boundary, we wish to expand the fixed point 
equations near the critical solution as we did in the dis- 
cussion of the undiluted model. Rather than working 
with equation |2| directly, it is more natural to use an 
equivalent infinite set of recursion relations for the mo- 
ments of P(h). These can be derived by multiplying both 
sides of equation ([2| by h n and integrating or by consid- 



h '(d) 



9-1 



tj tanh 1 ( tanh( J) tanh(/i^ 



(d-i). 



taken to the power n and averaged. Near P PM (h), 
expand this relation around small hf. 



9-1 

i=l 



1 



eiT(hi- -(1-T 2 )hi) + 



we 



(7) 



where T = tanh J and we have suppressed the depth 
superscripts. 

Sufficiently near the phase boundary, we expect the 
moments (h n ) to decrease exponentially with n and thus 
only a few leading order moments need be retained to 
extract the leading critical behavior at finite J. Taking 
powers of equation ([7]) and averaging, we find 

(h') = 2pT(h)- 2 -pT(\-T 2 )(h 3 ) 

(h 12 ) = 2pT 2 (h 2 ) + 2p 2 T 2 (h) 2 
(ti 3 ) = 2pT 3 (h 3 ) + 6p 2 T 3 (h) (h 2 ) 

to cubic order. We have specialized to the case q = 3 
in order to simplify the presentation; for q > 3 an ad- 
ditional term at cubic order is generated but the critical 
exponents remain the same. 

Near the phase boundary in the p — J plane, we can 
define a small parameter e by writing 2pT = 2- tanh J = 
1 + e. We will treat the fixed point equations to leading 
order in the e expansion. For e < the only real solution 
is paramagnetic: 



(h) = (h 2 ) = (h 3 ) = 0. 



(8) 



This solution is stable since its local Lyapunov exponents 
(e, e — In ^ , e — 2 In 4) are all negative. 

For e > this solution becomes unstable. We find 
two other ferromagnetically ordered solutions, which are 
linked by the symmetry h — > —h, and choose the positive 
one. This is 



(h) 
(h 2 ) 
(h 3 ) 



= 2 



1-T £l 



/2. 



T 



2 

— e 
T 



,3/2 



y/T(l -T)l + T 



(9) 



One can find the Lyapunov exponents of this stationary 
point analytically but the expressions are unenlighten- 
ing. We plot a typical case in Figure [2j Notice that the 
results (JsJ) are consistent with the assumption that (h n ) 
decreases exponentially with n. 

From Eq. Q one can read off the critical exponent 
(3 = 1/2 as the power of e in (h) ~ m. This is valid for all 
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FIG. 2: Lyapunov exponents of the ferromagnetic fixed point 
of the iteration equations for J = 1. Notice that for e < 0.13... 
they are all negative, signaling stability of the solution. At 
larger e, the stationary point becomes a focus before eventu- 
ally becoming unstable. 



T < 1 and sufficiently small e. The point T = 1 ( J = oo) 
is different and needs to be treated more carefully. As 
T — ► 1 the coefficient of e 1 / 2 in ^ vanishes, which implies 
that the J = oo critical exponent j3' > 1/2 while the 
divergence of the coefficient of e 3 / 2 means that (3' < 3/2. 
Indeed, from the exact solution of Section IV we will find 

p = i. 

For sufficiently large e > e c , the Lyapunov exponents 



become positive, signaling a loss of stability of the third 
order ferromagnetic solution (for the value J = 1, T = 
tanh(l) in Figure [2j e c = 0.137... ). This indicates that 
the first few moments flow to large scale and our trun- 
cation to cubic order fails. The value of e c decreases 
monotonically as T approaches 1 and to accurately find 
the fixed points we need to keep track of more moments 
of h in our iteration equations. At this point it is con- 
venient to switch to numerical solution of the full cavity 
equation ^ by population dynamics, as described in the 
Appendix. 

Having explored both above and below the critical 
point, we return briefly to the critical point at e = 0. 
Here, at linear order, there is a marginal flow near the 
paramagnetic fixed point. It is possible to analyze the 
truncated flow equations ^ at higher order to discover 
that the paramagnetic solution is indeed algebraically 
(rather than exponentially) stable, as one expects of a 
second order phase transition. With some additional 
algebra it is possible to carry a small applied field H 
through all of the above arguments at e = and show 
that the critical exponent (5 = 3 all along the p > p c 
phase boundary. 

The final important feature of the phase diagram is 
the presence of GM singularities throughout the p < 1, 
J > Jq region. That is, the density of LY zeros on 
the imaginary H axis of the partition function has an 
essential singularity like er a / lmH throughout this region 
due to the cumulative influence of rare large undiluted 



regions and there is therefore no gap. Equivalently, the 
real magnetization M ~ e~ a / H in a real applied field 
H. Although this can be seen from elementary rigorous 
arguments [BJ[7], it is difficult to detect either analytically 
or numerically at finite J. However in Sections |VI| and 
|VII| we will use the exact solution of the cavity equations 
at J = oo to exhibit these essential singularities explicitly 
and subject them to detailed study. 



IV. CLUSTER SERIES AT J = oo 

For the remainder of the paper, we will focus primarily 
on the J = oo part of the phase diagram of the model. 
We first review the classic argument due to Harris [T3] 
based on an expansion over connected clusters. This will 
lead to an exact series expansion for M(H) that we will 
independently rederive using the cavity approach in Sec- 
tion El 

Consider a cluster of n+ 1 spins connected by n bonds. 
For J ^> H ~ 1 (which is the meaning of the J = oo 
limit) each connected cluster behaves like a piece of fcr- 
romagnet. Indeed for H = there are two degenerate 
ground states, one with all spins pointing up and one 
with all spins pointing down. The first excited states 
are spin flips at energy ~ J above the ground states and 
their presence is negligible. Turning on a magnetic field 
H the degeneracy is broken and (if H is positive, say) the 
state with all spins pointing up is energetically preferred. 
Therefore the cluster will acquire a small magnetization: 



M n (H) = (n + 1) tanh((n + 1)H). 



(10) 



The total magnetization per spin is obtained by summing 
over all the clusters with their weights W n , corresponding 
to the number of clusters of size n per spin |19j 



M(H) 



n>0 



W n M n (H). 



(11) 



This equation has been studied before and results can 
be found in |13j for the magnetization, and |16j for the 
scaling law of the magnetization at the critical point. 
We will reproduce those results on the magnetization for 
completeness, but the main focus of this paper will be the 
density of Lee- Yang zeros and the solution of cavity field 
equations from which we will recover the known results. 

From the solution of the bond percolation problem on 
the Bethe lattice [12] the number per spin W n of clusters 
of bond size n is given by 



W n (p) 



. ((rc+l)(g-l))! 
■(n+l)\(n(q-2) + q)\ 



For simplicity we consider q — 3. Then p c = 1/2. We 
can easily obtain the asymptotic behavior of the magne- 
tization by using the asymptotics of W n as 



5/2 



,-nA(p) 



(12) 



where 



A(p) = In 



4p(l-p) 



(13) 



A(p) is the exponent governing the decay rate of the clus- 
ter sizes and it will appear often in the remainder of the 
paper. For p < p c — 1/2, A > and W n decreases expo- 
nentially. For p = p c = 1/2, A — and we have instead 
a power-law decay with exponent 5/2 (the exponent is 
independent of q) : 



5/2 



(14) 



This change in the asymptotic fall-off of the cluster dis- 
tribution at criticality is the reason for the change in the 
response to an applied external field at zero temperature. 

Indeed we can easily see how this works. For A > 
and small H we can write an asymptotic expansion: 



M(H) = W n (n+ l)tanh((n- 

n 

0(H 3 ) 



1)H) 



= H 



(1+P) 
1 - (q - l)p 



0(H 3 ) 



(15) 



which is linear in H. At the percolation threshold, how- 
ever, A = and (n 2 ) diverges. The expansion of tanh 
inside the first sum is unjustified. To find the first term 
in the asymptotic expansion of M (that we will derive 
in a formally correct way in Section VI) we use instead 



(121 



M{H) ~ Y, 



3/2 



tanh(ni7) 



(16) 



3 f°° 3 

— ^VH / dx x'Hanhx + O(H) 
2 V 7r Jo 



So the susceptibility diverges although there is no spon- 
taneous magnetization |13j . We now turn to a derivation 
of the above results from the cavity method. 



V. CAVITY APPROACH AT J = oo 

The cavity method for this system gives a probability 
distribution for the cavity fields which satisfies the fixed 
point equation (cf. Equation p|) 

P(h) = E e [ J| dh t P(hi)S (h~Y u ( h i + H ' Je i)) 
"* i=l V z=l ) 

(17) 

In the J — > oo limit, we can linearize the cavity biases u: 
u{h + H,Je l ) = {h + H)e l (18) 



At q = 3, the fixed-point equation ( 17 1 becomes 
P(h) = (l-p) 2 8(h) + 2p(l-p)P(h-H) 



P 2 J dh 2 P(h 2 - 2H)P(h-h 2 



This equation can be solved by defining the Laplace 
transform 



9(s) 



dhP(h)e 



- s h 



making sure to include the delta function at h — 0. The 
equation for g is quadratic 



= p 2 e- 2sH g( S ) 2 + (2p(l-p)e- sH -l)g(s) 
with solution 

„1Hs _ 2e tf»(p_p2)_ e 3l # S 



(19) 



(20) 

The second solution to ([19| is not physical. Even without 
inverting the Laplace transform all the properties of the 
solution can be extracted from g(s). For example the 
normalization condition, the zeroth moment, is 



dhP{h) = g(0) 



1 ifp<l/2; 
^ ifp>l/2. 



and the first moment is 
(h) = - 



dg 
ds 



2Hp 
l~2p 



(21) 



whose divergence at p = p c = 1/2 signals the ferromag- 
netic phase transition. For p > p c , P(h) loses normal- 
ization because a finite fraction of the cavity fields flow 
to infinity, just as in a percolating cluster distribution. 
Indeed, the divergent cavity fields are precisely those at- 
tached to spins in percolating clusters. Because these 
spins are connected to the positively biased boundary and 
the temperature is effectively zero, they spontaneously 
magnetize to M = 1 = tanh(oo). This provides the spon- 
taneous magnetization critical exponent: 

M(H = 0, J = oo,p) = 1 - Q—p- oc (p - p^ 1 

from which we read j3' = 1 in accord with the discussion 
following Eq. 

We now concentrate on the p < p c = 1/2 region at 
finite H. Consider the magnetization per spin 



M(H) = ^tanh(H + ^ u(hi + H, J a)) 



(22) 



i=l 



e,hi 



which is obtained by averaging over the disorder and the 
distribution of h. It is straightforward to show that for 
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small h and H we obtain the results of Equation ( 15 I 
Indeed: 



The a n are given by the series expansion of the square 
root in (p0|: 



M{H) = iH + J2^i(hi+H)) +0(H 3 ) (23) 

\ »=1 / e,h 

= (l+pq)H+qp{h) + 0(H 3 ); (24) 



(25) 



Now substitute ( plj ) and simplify 



in accordance with (15 1 



We now reconstruct the full probability distribution 
P(h) exactly. We expand the function g(s) as a series in 



(26) 



n>0 



which defines the coefficients a n . P(h) is now given by 
the inverse Laplace transform, 



P(h) = J2^ n S(h-nH). 



(27) 



n>0 



3=0 



{Ap(l~p)) n + 2 (-l)" +1 r(3/2) 



2p 2 



r(n + 3)r(-?i - 1/2) 



(28) 



That is, P(h) is a comb of delta functions at integer mul- 
tiples of H with a decaying envelope. The large n behav- 
ior of the envelope is 



a n = 4 ( 1 -P) 2 n -3/ 2e -A(p)n 
V7T 



(29) 



where A(p) is defined in (13 1. It is not surprising that 



the same asymptotics governs both P(h) and W n . 

Finally, to connect directly with the previous section 
let us compute the exact magnetization of a spin as a 
function of applied field H. Evaluating the cavity magne- 
tization equation ( 22 I using the cavity field distribution 
p7), we find 



M(H) = ^n\p i (l-p) q - i «mi-"«m i tanh((j + l + mi + ..- + m i )fl) 



(30) 



which naturally expands as a series in tanh nH. At q — 3, 
we can evaluate all the coefficients in this series to find 
that indeed they are identical to the coefficients nW n -i 
of Eq. (111. Thus, at J = oo the cluster series and the 



cavity method produce identical results for the magneti- 
zation. 

Having established the equivalence of the two solu- 
tions, we now return to the analysis of the series for the 
magnetization. In the following two sections we will ex- 
tract the critical behavior near p = p c and, by analytic 
continuation to imaginary H, the GM singularity in the 
density of LY zeros. 



VI. 



AN INTEGRAL REPRESENTATION FOR 
THE MAGNETIZATION 



Despite its simplicity, the expansion of Eq. ( 1 1 



is an 

exact result for the magnetization which can be analyti- 
cally continued to imaginary values of the magnetic field. 
However, the representation of M(H) as a sum in (11) 
is not best suited for this purpose. An integral repre- 
sentation would be preferable. To obtain it we write the 



Laplace transform of the function tanh x 



M = 



dxe sx tanh : 



1> 



2 + s 
4 



(31) 



where ip is the digamma function. The function /(s) has 
simples poles only at the negative even integers and thus 
we can invert the transform and write 



tanh x 



ds 
2tH 



e sx f{s), 



(32) 



where B is any Bromwich path lying to the right of all 
poles of f(s), that is to the right of the negative real 
axis. Inserting into (111, we can invert sum and integral, 
provided 

|4p(l -p)e siI \ < 1. (33) 
The resulting expression, valid for |argiJ| < 7r/2, 

M(H) = 3(1 -p) 3 / — /(a) 
Jb 2m 



V 2(n+1)! e s{n+1 ^ H (v(l - v)Y 
^-(n + 3)!n! [P[ P)) 



n>0 



7 



can be written in closed form by performing the sum. erating function of the probability W n . For the Bcthc 
This amounts to calculating the derivative of the gen- lattice the generating function is [20] 



<t>{x) = ^2w n x n 

n>0 

= -2(l-p) 3 a; - 3 (8(l-yW) + 4(2yr^-3)e + 3C 2 ) (34) 

where £ = 4p(l — p)x has been defined for convenience. By means of this function we can perform the sum inside the 
integral to obtain 

3(1 -pf f ds „^_ 2aH 

^ (35) 



x f(j)(l -p)e sH 4p(l -p)e sH + 1 - 3p(l - p)e sH ^j 



We simplify this expression by considering the analytic 
structure of the integrand. The function /(s) has simple 
poles only at the negative even integers (0,-2,-4,...), while 
the rational expression has a series of square root cuts at 
■ s n = ^tt~ *27m/ H. We close the contour with a semi- 
circle at infinity on the right (for Re s > 0) on the first 



Riemann sheet. We then deform the contour to coincide 
with the edges of the cuts. At this point only the discon- 
tinuity across the cuts contributes to the final result. For 
aesthetic reasons we finally shift the value of s by ^jip-, 
the real part of the origins of the cuts. 
The resulting expression is 



M(H) 



y, 8(1 -p)' 
n=—oc 



(36) 



At this point it seems we have traded a sum of func- expanding the integrand in powers of H we find 



tions (16 1 with a series of integrals that we cannot eval- 
uate. This looks like a step backward in the quest for 
a useful result! However, after thinking about the pro- 
cedure we have performed, we recognize that this is a 
Poisson summation-like duality on the original equation 
(111. The terms in the sum are higher and higher har- 



monics of the result (this is particularly evident, as we 
will see shortly, for imaginary H). 



The scries in n in (36 1 is dual to the series in (111 so 



that when the first converges rapidly the second does not 
and vice versa (for H on the real axis). In the interesting 
regime, close to the percolation threshold (111 converges 
slowly and the first term (n = 0) of (36 1 gives the leading 
term in the expansion in (p — p c ) and H — > 0. 



Let us now see how we can recover Eq. ( 16 1 in a clean 
way. At the critical point p = 1/2, we have A — and 
so s* = i2irn/H. For H — > all the cuts except that 
corresponding to n = go to infinity and we can keep 
only the n = term in the series ( 36 1 . Moreover, by 



M(H) ~ -VH I Vsf(s) + 0(H) 

7T 



" 3/2 tanh x + O (H) 



which coincides with Eq. (16 1 



To recapitulate, the magnetization is given by an inte- 
gral of the discontinuous part of the generating function 
4> of the cluster distribution W n with the Laplace trans- 
form of the function tanhx. For the Bethe lattice the 
generating function can be written down explicitly and 
the calculations can be carried to the end. In the per- 
colation limit the cut on the real axis gives the greatest 
contribution to the sum. 



VII. DENSITY OF LEE- YANG ZEROS AT J = oo 

In this Section we will find the density of Lee- Yang 
zeros p at J = oo. These are the zeros of the partition 
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FIG. 3: The density of Lee- Yang zeros as a function of the 
imaginary field 6 at p = f/4. In red the smoothed p am , in 
orange the sum of the first three harmonics in ( 36 1 (terms 
n — ±I,±2,±3) and in blue 



-p sm - The figure suggests (in 
agreement with the discussion in the text) that the sum of all 
the harmonics (with n 7^ ) builds a sum of delta functions 
Eq. (371 minus p am in (40 1. 



function as a function of the external magnetic field H, 
for imaginary H — id. Instead of solving the equation 
Z(H) = directly, we rely on the relation 

p(6) = -ReM(i9 + 0+). 



To get an idea of how this function looks we recall 



Re tanh(i6» + 0+) = it ^ 6(6 - |(2m + 1)). 



m— — oo 



From this wc find 



p(°) = EE^»( n+1 M n+1 ) < '-?( 2m + 1 )) 



m n>0 



2m + 1 
2n + 2 



(37) 



so the zeros are located at all the - 2£ ^- rational multiples 

even r 

of it, with multiplicities given by the WnS [21]. This 
is a singular distribution with an accumulation point at 
6 = 0: our task is now to smooth it by using the Poisson- 
dual integral representation (36 1 obtained in the previous 
section. 

The expansion over the cuts, for imaginary H becomes 
an expansion in higher and higher harmonics (see Figure 
|3| of p. Selecting the term with n — in ( 36 1 , gives the 
function smoothed to the lowest degree: 



ft m («) = f <We se - 1 (l - \e se ^j e- 2se Im f{-is - iA(p)/8 + 0+). (38) 

This expression simplifies since from the definition of / 



Im f (— iz + + ) — I dx sin(za;) tanh x = — . (39) 

' 2sinh7rz/2 



So we find the smoothed density of LY zeros 



4(1 - P f 



Psm(0) = — — — / dsVe se -1(1- - t e sB ) ^- — — . (40) 

irp Jo V 4 J smh(f(s + A/0)) 



e 



The different profiles for this density can be seen in 
Figure [4] Here the GM phenomenon is evident: even at 
p < p c = 1/2, p S m(0) is strictly positive for any non- 
zero 9; there is no gap in the distribution. This effect is 
due to the presence of rare large clusters. The asymp- 
totic expansion of the density at small 8 can be found by 
expanding the integrand 



which is p = p c , where it shows the critical square root 
cusp in the magnetization. However, let us remark that 
Eq. (40) is the smoothed part (in the sense of distribu- 



tionsjfbr all values of 9 and not only for small 9. 



p{e) . *W\-r)V ee-w. 



(41) 



This expansion is uniformly valid at the point A = 0, 



Let us now make a few qualitative remarks on the 
asymptotic expansions for M(H) and p which apply in 
principle to all lattices. From the integral representa- 
tion we observe that there is no Stokes phenomenon for 
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FIG. 4: The density of Lee- Yang zeros as a function of the 
imaginary field 9 at p — 1/10 (above) and at criticality p = 
1/2 (below). 



VIII. SUMMARY 

The diluted Ising ferromagnet on a Bethe lattice is a 
tractable model that beautifully illustrates many of the 
key physical features of short-ranged disordered systems. 
In this paper, we have attempted to present a unified 
analysis of the model in the framework of the cavity 
method, from which we derive both well-known elemen- 
tary results about its phases and non-trivial features such 
as GM singularities and the infinite coupling critical ex- 
ponents. 

In particular, the ferromagnetic phase boundary lies 
in the mean-field universality class (6 — 3) at any di- 
lution above the percolation threshold. At this thresh- 
old however, the ferromagnetic critical coupling diverges 
(J — > oo ) and our closed form solutions for the cavity dis- 
tributions in this limit reveal that the critical behavior 
is governed by the percolation of the underlying lattice 
(6 = 2). Linear stability analysis of the flow of the cavity 
moments near criticality naturally reveals the Lyapunov 
exponents and the associated correlation depth of the 
stable phases. 

Furthermore, at infinite coupling we have explicitly ex- 
hibited the essential Griffiths-McCoy singularities in the 
magnetization for all p < p c , where there is no sponta- 
neous magnetization. By an harmonic rcsummation of 
the exact magnetization, we found the smoothed density 
of LY zeros exactly and conjectured its relation to the 
real part of an appropriate terminant of the asymptotic 
series for the magnetization. 



KeH > 0, i.e. the asymptotic approximation for M(H), 



M L (H) = a k H 2k+1 ((n + l) 2k+2 )+0 (H 



2L+l\ 



(42) 



fc=0 



(where a^'s are the coefficients in the series expansion 
for tanhx) is valid for all |argi/| < tt/2. Naively, sub- 
stituting H = i9 + + into this expansion, we obtain a 
purely imaginary result for any L and we might spec- 
ulate that p = 0. However, the expansion (42 1 is only 
asymptotic, since ((n+ l) fe ) ~ k\e~ kA . This means that 
we cannot take L — > oo but rather truncate the series 
at the L ~ H/A where the remainder is smallest. The 
remainder is never actually zero but it is exponentially 
small in l/|i?|. A good quantitave approximation can be 
found by using the "terminant" |17j of the asymptotic 
expansion. The terminant is indeed cx e~% A /' H ' and it is 
not purely imaginary for H = i6 + + . Thus, as a gen- 
oral rule we expect the real part of the terminant of the 
asymptotic expansion of M(H) represents the density of 
LY zeros in the subcritical region. 
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APPENDIX: NUMERICAL METHODS 

In the cavity framework, all of the statistical observ- 
ables of a model can be derived from the cavity field 
distribution P(h). This distribution is the fixed point of 
the flow of the cavity equation ([2]). While we can solve 
this equation analytically in certain limits, we rely on 
a simple iterative numerical algorithm called population 
dynamics for many of the finite coupling results. Popula- 
tion dynamics and its more sophisticated variants appear 
in, for example, |18j . 
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FIG. 5: Cavity field distribution at p — 0.5, J — oo critical 
point with H = 1. This log-log plot shows agreement with the 
asymptotic form P(h) ~ h~ 3 ^ 2 found in the exact solution. 



The algorithm works as follows: we represent the dis- 
tribution P(h) by a finite population of N pop fields hi. 
This population is initialized from an appropriate uni- 
form distribution and then iterated as follows: 

1. Select q—l fields hi randomly from the population 
and q — l random e.j . 

2. Use ^ to calculate the cavity field ho on a spin 
sitting below the q—l spins selected above. 

3. Randomly replace one element of the population 
with ho. 

4. Repeat until convergence in some measure of the 
population, for example the cavity magnetization 
(tanh(/i)). 

In practice, this procedure converges quickly deep in ci- 
ther the ferromagnetic or paramagnetic phase but slows 
near the phase transition. We illustrate some typical re- 
sults below for q — 3. 

Even at the percolating critical point, when the ex- 
pected cavity distribution develops a long tail and diver- 
gent moments, this procedure works. Figure [5] shows the 
numerically determined cavity field distribution for the 
p = 0.5, J = oo critical point with applied field H = 1. 
As noted in Section [V] the exact solution is a comb of 
delta functions at h — nH, n € N with weights decay- 
ing asymptotically as a power law a n ~ n~ 3 / 2 . The 
numerical solution concentrates on integer fields with a 
power-law tail consistent with the exponent —1.5. 

In general, the form of the cavity field distribution at 
finite p and J is only obtainable numerically. Figure [6] 
shows the numerically determined distribution at small 



applied field H on three different points in the p — J plane 
near the p = 0.75, J = 0.80 critical point. These distri- 
butions are typical and illustrate the dramatic increase 
in susceptibility on the ferromagnetic side of the phase 
boundary. 
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FIG. 6: Cavity field distributions as function of p near the 
p = 0.75, J = 0.80 critical point at small applied field H ~ 
* The blue curve is in the paramagnetic regime, the green 



10" 



is in the critical phase and the red in the ferromagnetic one. 
Notice the dramatic increase in the response of the cavity field 
distribution to the applied field on the ferromagnetic side of 
the phase boundary. 
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FIG. 7: Critical magnetization in an external field along the 
J — <x line. The slope of the log-log curves indicates the 
critical exponent 8 associated with each of the four points 
p — 0.25, 0.5, 0.75, 1 along the phase boundary. We find 5 = 1 
for p = 0.25 < p c = 0.5, 8 = 1/2 for p = 0.5 and 8 = 1/3 for 
p = 0.75 and 1. 



Finally, Figure [7] confirms numerically the critical ex- 
ponents derived using the moment flow analysis of Sec- 
tion Ull] and the exact solution of Section II VI 
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